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Based on an ensemble Monte Carlo analysis, we show that Coulomb interactions play a dominant role in bound- 
to-continuum terahertz quantum cascade lasers and thus require careful modeling. Coulomb interactions enter our 
simulation in the form of space charge effects as well as Coulomb scattering events. By comparison to a full many- 
subband Coulomb screening model, we show that simplified approaches produce considerable deviations for such 
structures. Also the spin dependence of electron-electron scattering has to be adequately considered. Moreover, we 
demonstrate that iterative Schrodinger-Poisson and carrier transport simulations are necessary to correctly account for 
space charge effects. 

PACS numbers: 42.55.Px, 73.63.Hs, 78.70.Gq 



I. INTRODUCTION 

For terahertz quantum cascade lasers (QCLs), two types of 
structures play a major role: the resonant-phonon (RP) and 
the bound-to-continuum (BTC) design. In RP structures, ef- 
ficient depletion of the lower laser level is achieved by tun- 
ing the corresponding transition to the longitudinal-optical 
(LO) phonon energy (36 meV in GaAs). By contrast, BTC 
designs are based on minibands, enhancing the influence of 
Coulomb interactions in two ways. First, the close energetic 
spacing of the miniband levels favors electron-electron (e-e) 
over LO phonon scattering in the carrier transport. Also, the 
large spatial extent of the minibands across many wells, to- 
gether with the localization of the positively charged donors 
in typically a single well, leads to considerable conduction 
band bending due to space charge effects.' Thus, a careful 
modeling of Coulomb interactions is necessary to analyze car- 
rier transport in BTC QCLs, and here specifically the role of 
Coulomb interactions in such structures. However, e-e scat- 
tering is much more computationally demanding than single- 
electron processes like electron-phonon interactions, hamper- 
ing its inclusion in quantum mechanical simulations of QCLs 
beyond the mean-field approximation.^ The numerical load is 
further increased by the large spatial extent of the minibands. 
Due to its efficiency, the semiclassical ensemble Monte-Carlo 
(EMC) method is well-suited for investigating BTC struc- 
tures. However, EMC simulations have up to now typically 
focused on terahertz RP structures, ^"^ while only few results 
are available for equivalent BTC or related chirped superlat- 
tice designsi^iii 

We present an EMC simulation tool optimized for the sim- 
ulation of terahertz BTC structures. We employ a full many- 
subband (MS) screening model for electron-electron scat- 
tering and include the exchange effect for the parallel-spin 
colUsions. Space charge effects are adequately considered 
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by performing iterative Schrodinger-Poisson (SP) and carrier 
transport simulations, yielding so-called self-self-consistent 
solutions.' This tool allows us to properly assess the role of 
Coulomb interactions for the carrier transport and gain in a 
typical terahertz BTC design. Widely-used approximations, 
such as non-iterative simulations, using simplified screening 
models, or neglecting the spin dependence of e-e scattering, 
usually work well for RP structures, but can lead to significant 
deviations for terahertz BTC QCLs, as shown in this paper. 



II. METHOD 

The simulation tool, consisting of a three-dimensional 
EMC and an SP solver, allows for a self-consistent analysis of 
the carrier transport and optical gain in the QCL structure.'— 
The SP solver yields the subband eigenenergies and wave 
functions, needed as an input for the semiclassical EMC car- 
rier transport simulation. The obtained electron distribution 
in the structure gives rise to space charge effects, resulting in 
conduction band bending and altered subbands. Thus, itera- 
tive runs of SP and EMC simulations are necessary to obtain 
convergence, corresponding to self-self-consistent solutionsii 

All essential scattering mechanisms are accounted for 
in the carrier transport simulation. Included is elastic 
electron-impurity (e-i) and interface roughness scattering, as 
well as inelastic interactions of electrons with acoustic and 
longitudinal-optical (LO) phonons. Nonequilibrium phonon 
effects are also taken into account. '^ Being evaluated as a two- 
electron process, e-e scattering assumes a special role in the 
EMC simulation and has the highest computational complex- 
ity. Its implementation is more closely discussed in Section 
IIIBI For all scattering processes, Pauli's exclusion principle 
is taken into account.'^ Periodic boundary conditions are em- 
ployed, i.e., electrons leaving the device on one side are auto- 
matically injected into the equivalent subband on the opposite 
side.-' 

Coulomb interactions enter our simulation in the form of 
SDace charge effects as well as individual e-e and e-i scatter- 
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ing events. As discussed above, space charge effects and e-e 
scattering play an especially important role in THz BTC struc- 
tures due to the formation of minibands. Thus, the implemen- 
tation of these mechanisms will be discussed in detail in the 
following. 



A. Space charge effects 

The subband wave functions \\t„ (z) and eigenenergies E„ 
for the simulated QCL structure are obtained by solving the 
SP system^i 
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Here, a position independent permittivity e is assumed. Fur- 
thermore, m* (z) and (z) are the electron effective mass and 
doping concentration of the structure, e and h are the elemen- 
tary charge and the reduced Planck constant, and n2D,« is the 
electron sheet density of level n. The self-consistent poten- 
tial is given by V (z) — Vo{z) — ecp(z), where Vb is the con- 
duction band profile and (p is the electric potential due to the 
space charge profile. For QCLs, the subband energies and 
wave functions are typically computed for a single central pe- 
riod of the structure; the subbands in other periods are then 
obtained by appropriate shifts of the solutions in energy and 
position. For the first run of the SP solver, a thermal occu- 
pation of the subbands according to Fermi-Dirac statistics is 
assumed.''* To obtain self-self-consistent solutions, the sub- 
band occupations are subsequently extracted from the EMC 
analysis, which is carried out alternately with the SP simula- 
tion until mutual convergence is obtained. 



B. Electron-electron scattering 

In the EMC simulation, e-e scattering is implemented as a 
two-electron processi '^''^ An electron in an initial state jik), 
i.e., subband / and in-plane wave vector k, scatters to a final 
state |y'k'), accompanied by a transition of a second electron 
from a state \iok.o) to |yoko). The total scattering rate from 
ik) to a subband j is then obtained by the Fermi golden rule. 
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with the electron effective mass m*, cross section area A and 
carrier distribution function /,(, in subband /q. 9 is the angle 
between g = ko — k and g' = kg — k', and Q = k — k' (with 
G = |Q|) denotes the exchanged wavevector. 

Different approaches with varying degrees of complexity 
exist to compute the transition matrix element Mugjjg from 



the bare Coulomb matrix elements, 
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First, the screened Coulomb matrix elements V,; (Q) are 

obtained from y^^jj^^ {Q) by applying a more or less sophis- 
ticated screening model. In the random phase approximation 
(RPA), they are found by solving the equation systen>ii 
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Here, n,„„ {Q) is the polarizability tensor, given in the long 
wavelength limit {Q 0) by 
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For collisions of electrons with parallel spin, interference 
occurs between V,? ,-, and the 'exchange' matrix element 
y? . .J^ Accounting for this exchange effect, the magnitude 
squared of the transition matrix element Mugjjg is then given 

b\"'" 
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where = pp = 1/2 are the probabilities for antiparallel and 
parallel spin collisions, respectively, and 



2g^+gl±2g(g^ 



-glf" cose 



1/2 



(7) 



with g = |g| and ^2 ^ 4m* (£, + - Ej - /fil 

Commonly, simplified screening models are used to avoid 
the numerical load associated with solving Eq. (01).^'^ Further- 
more, often the exchange effect is neglected when calculating 
Mii^jj^^}^ As discussed in SectionU e-e scattering plays a par- 
ticularly significant role in BTC structures. Thus, the valid- 
ity of such approximations in BTC designs deserves special 
scrutiny, as shown below. 



III. RESULTS AND DISCUSSION 

We present simulation results for a 3.5 THz BTC design,— 
demonstrating the importance of Coulomb interactions on the 
carrier transport and material gain in such structures. In par- 
ticular, we investigate the validity of neglecting the exchange 
effect and employing simplified screening models for the eval- 
uation of e-e scattering. Also the importance of self-self- 
consistent simulations for a proper inclusion of space charge 
effects is discussed. In our EMC simulation, all scattering 
mechanisms are evaluated self-consistently. Only interface 
roughness (IR), whose parameters are hard to measure and 



3 



depend critically on the growth conditions, has to be de- 
scribed in terms of a phenomenological model. For the 
IR mean height and correlation length, we use typical values 
of A = 0.12nm, r= 10nm,ii*^ 



A. Space charge effects 

In Fig. [T] the simulated conduction band profile and spec- 
tral gain of the investigated BTC QCL is shown for a lattice 
temperature IL = lOK at the design bias of 2.5kV/cm, where 
the maximum gain is obtained in the simulation. Solid lines 
indicate fully self-self-consistent results; here, the SP and 
EMC simulations are carried out iteratively until convergence 
is obtained. Dashed lines indicate the results obtained by solv- 
ing the SP system once in the beginning assuming a thermal 
carrier distribution,'^ and then performing a self-consistent 
EMC simulation, which is a quite common approach . 
For comparison, also the conduction band profile and spec- 
tral gain is shown as obtained with deactivated Poisson solver 
(dotted lines), i.e., when no space charge effects are included. 
In all cases, e-e scattering is evaluated taking into account the 
exchange effect as well as screening, which is considered in 
RPA by repeatedly solving Eq. ^ to account for changes in 
the carrier distribution during the EMC simulation. 

The results shown in Fig. [T] indicate that the proper in- 
clusion of space charge effects greatly affects the simulation 
outcome for the investigated BTC structure. A significant 
conduction band bending is observed for activated Poisson 
solver (solid and dashed lines in Fig. \lla)), changing the sub- 
band eigenenergies and wave functions. Comparison to ex- 
perimental data shows that only the simulation results with 
space charge effects included are in line with experiment. The 
gain profile obtained with deactivated Poisson solver is ex- 
ceedingly broad, flat and low, see Fig. [Ttb) (dotted curve). 
The low peak gain of around 9cm^' at around 4.2 THz is in 
sharp contrast to the experimentally observed lasing of the 
structure at 3.5 THz. Also the large gain bandwidth is not 
in accord with electroluminescence measurements, yielding 
full width at half maximum (FWHM) widths of clearly below 
1 THz.'^ The inclusion of space charge effects leads to a real- 
istic gain profile centered around 3.5 THz in agreement with 
experiment. The FWHM widths of 0.65 THz (solid curve) and 
0.66 THz (dashed curve) agree reasonably well with the ex- 
perimental value of 0.85 THz at lOK, extracted from electro- 
luminescence measurements.'^ Still, the two gain profiles are 
somewhat different, with peak gain values of 24.2 cm^' for 
the self-self-consistent and 21.2 cm^ ' for the non-iterative ap- 
proach. This illustrates the importance of self-self-consistent 
simulations for terahertz BTC designs, where space charge ef- 
fects tend to play a more pronounced role than in equivalent 
RP structures. 



B. Exchange effect 

There are two common approaches to implement e-e scat- 
tering without explicitly considering the spin dependence. 
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FIG. 1. (Color online) (a) Conduction band profile and probability 
densities, as obtained for thermally occupied subbands (dashed lines) 
and by a self-self-consistent simulation (solid lines). For compari- 
son, also the conduction band profile without space charge effects 
included is displayed (dotted line), (b) Simulated spectral gain ver- 
sus frequency for the cases shown above. 
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FIG. 2. (Color online) Simulated spectral gain vs frequency, as ob- 
tained by fully taking into account (solid curves) and ignoring the 
exchange effect (dotted curves), and by ignoring parallel spin colli- 
sions (dashed curves). For comparison, the result obtained with no 
e-e scattering is also shown (dash-dotted curves), (a) Tl = lOK; (b) 
Tl = 90 K. 



One method, which tends to overestimate the exchange ef- 
fect, is to completely neglect the parallel spin collisions , '^''^ 
implying pa = 1/2, = in Eq. Another common 

approach is to ignore the spin dependence. Parallel spin 
collisions are then treated the same way as antiparallel spin 
contributions,'** which corresponds to p,, = 1, pp = in Eq. 
(|6]l. In moderately doped RP structures at elevated tempera- 
tures, where the carrier transport is dominated by LO phonon 
scattering, the contribution of the exchange effect is usually 
negligible. In BTC designs, based on minibands with closely 
spaced energy levels, e-e scattering plays a more pronounced 
role." 

Figure |2] contains the self-self-consistently simulated gain 
spectra at 7l = lOK and TL = 90 K, taking into account 
screening in the RPA. Results are shown for pa = Pp = 1/2 
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FIG. 3. (Color online) Simulation results for (a) relative subband 
occupations and (b) subband temperatures, as obtained by fully tak- 
ing into account (x-marks, solid curves) and ignoring the exchange 
effect (plus signs, dotted curves), and by ignoring parallel spin colli- 
sions (circles, dashed curves). The lines are guide to the eye. 



FIG. 4. (Color online) Simulation results for electron dwell times, 
as obtained by fully taking into account (x-marks, solid curves) and 
ignoring the exchange effect (plus signs, dotted curves), and by ig- 
noring parallel spin collisions (circles, dashed curves). The lines are 
guide to the eye. 



(solid curve), /?« = 1, = (dotted curve), /?„ = 1/2, /?p = 
(dashed curve), and Pa= Pp = (dash-dotted curve). The last 
case, which corresponds to completely neglecting e-e scatter- 
ing in the simulation, yields two narrow gain spikes around 
2.8 and3.6THz, largely deviating from the experimental elec- 
troluminescence measurements.'^ This illustrates the impor- 
tance of e-e scattering for such structures. As discussed in 
Section lTlI Al for Tl = lOK the full simulation with exchange 
effect included yields a realistic gain profile, see solid curve in 
Fig. |2a). Ignoring exchange (dotted curve) leads to an over- 
estimation of the scattering, resulting in a peak gain reduction 
by 8%, and an increase of the FWHM gain width by 11%. On 
the other hand, completely neglecting parallel spin collisions 
(dashed curve) leads to a gain width reduction of 24%, and 
an increase of the peak gain by 15%. For TL = 90 K, see Fig. 
I2b), the gain gets somewhat broadened and lowered in agree- 
ment with experiment;" '^ also here, ignoring the exchange 
effect or parallel spin collisions has similar effects on the sim- 
ulated gain profile as discussed for TL = lOK. 

In Fig. |3] the relative occupations and the electron tempera- 
tures are shown at Tl = lOK for the eight energy levels within 
a miniband, characterized by their corresponding eigenener- 
gies (compare Fig. [TJa))- The electron temperatures are ex- 
tracted from the (generally non-thermal) carrier distributions 
in each subband by a least square fit. While the subband oc- 
cupations do not depend much on the implementation of the 
exchange effect, the electron temperatures show a stronger de- 
pendence. For example, the extracted temperature for the 8th 
level (at 39.3 me V) varies from 121. 6K to 146. 7K, depend- 
ing on the implementation of the spin dependence. The rela- 
tive insensitivity of the population and the strong dependence 
of the electron temperature and gain on the implementation 
can be understood by looking at the average dwell time of 
an electron in a given subband, which is the inverse of the 
outscattering rate from this subband. In Fig. H] the dwell 
time is shown for the eight energy levels, again characterized 
by their eigenenergies. The scattering is lowest, i.e., the dwell 
time is highest, when parallel spin collisions are ignored in the 
simulation (dashed curve). This is consistent with the reduced 
bandwidth and enhanced peak value of the corresponding gain 



profile in Fig. |2|a), which are directly related to the outscat- 
tering rate.— On the other hand, the electron dwell time in- 
creases by a similar factor for all the subbands, thus explaining 
the relative insensitivity of subband occupations to the chosen 
implementation of Eq. 

C. Screening 

Due to the computational effort involved in the RPA, 
screening is commonly taken into account using simplified 
models rather than solving Eq. directly. The screen- 
ing can for instance be considered by introducing a screen- 
ing wavenumber in Eq. (|3]l, i.e., replacing the prefactor 
e^/ {2eQ) by e^/ [2e{Q + qs)]. In single subband models, 
is obtained from Eq. dUi by assuming that screening is caused 
only by a single subband, e.g., the ground state. ^-^ A modified 
approach, which has been shown to yield improved results for 
RP structures, is the modified single subband model^ with 

where / sums over the subbands in one period. 

In Fig. |5] the intrasubband Coulomb matrix elements Vmi 
(Fig- Ha)) and Vnn (Fig. |3b)), as well as the intersub- 
band element Vmi (Fig. |5|c)) are shown as a function of the 
wavenumber Q. Here, 1 and 2 denote the upper laser level at 
22.2meV and the level directly above at 23.4meV, see Fig. 
[Tta). Displayed are the screened matrix elements based on 
the RPA (solid lines) and the simplified model according to 
Eq. (Is) with qs = 0.0237 nm^' (dashed lines), as well as the 
bare matrix elements defined in Eq. ^ (dotted lines). As can 
be seen from Fig. EJc), in the simplified screening model the 
intersubband elements approach zero for small wavenumbers, 
in contrast to the exact implementation of the RPA. A better 
approach consists in applying the simplified screening model 
only to the intrasubband matrix elements, and to treat inter- 
subband scattering as unscreened. 

Figure |6] contains the self-self-consistently simulated gain 
spectra at 7l = lOK and TL — 90 K. In contrast to Fig. |2] 
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FIG. 5. (Color online) Screened and unscreened Coulomb matrix 
elements for the structure shown in Fig. [TJa). (a) Vmi ; (b) V1212; (c) 
Vll22- 
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FIG. 6. (Color online) Simulated spectral gain vs frequency, as ob- 
tained by taking into account screening in the RPA (solid curves), 
and using the modified single subband model for all matrix ele- 
ments (dashed curves) or for the intrasubband elements only (dot- 
ted curves), i.e., treating intersubband elements as unscreened, (a) 
Jl = lOK; (b) Jl = 90 K. 



the exchange effect is included, but different screening mod- 
els are used. The reference curve based on the exact evalu- 
ation of the RPA (solid curves) agrees well with experiment, 
see Section IIII A I Applying the simplified screening model 
to all matrix elements overestimates the screening of the in- 
tersubband elements, compare Fig. |5jc), and thus results in 
an underestimation of scattering. The resulting spectral gain 
profile at lOK (dashed curve in Fig. |6ta)) features a 25% en- 
hanced gain peak and an excessively narrow FWHM width of 
0.43 THz, as compared to an experimental value of 0.85 THz. 
On the other hand, completely ignoring the screening effect 
for the intersubband matrix elements overestimates the inter- 
subband scattering, thus resulting in a lowered and broadened 
gain profile (dotted curves). As can be seen from Fig. |6jb), 
the simulation results at Tl = 90K are affected in a similar 
way by using above discussed approximations. In Fig. [T] the 
obtained relative subband occupations and electron tempera- 
tures are compared at TL = lOK for the different implemen- 
tations of screening. Although the simulated gain in Fig. |6] 
greatly depends on the applied screening model, the occupa- 
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FIG. 7. (Color online) Simulation results for (a) relative subband 
occupations and (b) subband temperatures, as obtained by taking 
into account screening in the RPA (x-marks, solid curves), and using 
the modified single subband model for all matrix elements (circles, 
dashed curves) or for the intrasubband elements only (plus signs, dot- 
ted curves). The lines are guide to the eye. 



tions of the miniband levels in Fig. [TJa) are quite insensitive 
to the chosen implementation, similarly as in Fig. [3] This is 
again consistent with the fact that the average dwell time of 
an electron in a level, shown in Fig. [8j strongly depends on 
the chosen screening model, but is changed by a similar factor 
for all subbands. We attribute this to the quasi-continuous na- 
ture of the minibands in the BTC structure. Contrariwise, for 
RP QCLs, where the energetic separation of the levels greatly 
varies, simple approaches like the single subband screening 
model have been shown to considerably affect the simulation 
results for the level occupations.^-^ Here, only a single sub- 
band, typically the ground level, is considered.? However, it 
should be mentioned that also for the investigated BTC struc- 
ture, this approach would yield somewhat less accurate results 
than evaluating Eq. as done in our simulations. The ki- 
netic electron distribution within each subband, represented 
by fitted electron temperatures in Fig. Illb), shows a mod- 
erate dependence on the screening model. For example, the 
extracted temperatures in the 8th level (at 39.3 me V) ranges 
between 127.9 K and 1 39.0 K for the different screening mod- 
els. Here again, more significant deviations from the RPA 
result would be obtained by using a single subband screening 
model rather than Eq. ([8]|. 



IV. CONCLUSIONS 

Employing an EMC analysis, we have investigated the 
role of Coulomb interactions in THz BTC structures, and as- 
sessed the validity of common approximations for modeling 
Coulomb effects. We have shown that space charge effects 
lead to considerable conduction band bending, having a sig- 
nificant influence on the obtained wave functions and the en- 
ergy level alignment. These mechanisms are properly ac- 
counted for in a self-self-consistent approach, i.e., by itera- 
tively performing SP and EMC simulations until mutual con- 
vergence is achieved. Ignoring these effects leads to exces- 
sively lowered and broadened gain profiles, in sharp contra- 
diction to experimental results. Also e-e scattering between 
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FIG. 8. (Color online) Simulation results for electron dwell times, as 
obtained by taking into account screening in the RPA (x-marks, solid 
curves), and using the modified single subband model for all matrix 
elements (circles, dashed curves) or for the intrasubband elements 
only (plus signs, dotted curves). The lines are guide to the eye. 



the closely spaced miniband levels plays an important role 
for the carrier transport and the spectral gain profile. Com- 
mon approximations in implementing e-e scattering, like ne- 
glecting the spin-related exchange effect or using simplified 
screening models, results in considerable deviations for the 
obtained gain profile. Even refined approximations, like the 
modified single subband screening model, introduce an error 
of approximately 25% for the peak gain in the investigated 
structure, as compared to simulations based on the full RPA. 
Combined with ignoring the spin dependence, the deviation 
can easily exceed 30%. With space charge effects and e-e 
scattering properly accounted for, the EMC analysis is shown 
to yield meaningful results for the investigated BTC design, 
which are found to be in good agreement with experiment. 
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